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ABSTRACT 

We present results of a population synthesis of millisecond pulsars from the 
Galactic disk. Excluding globular clusters, we model the spatial distribution of 
millisecond pulsars by assuming their birth in the Galactic disk with a random 
kick velocity and evolve them to the present within the Galactic potential. We 
assume that normal and millisecond pulsars are standard candles described with 
a common radio luminosity model that invokes a new relationship between ra- 
dio core and cone emission suggested by recent studies. In modeling the radio 
emission beams, we explore the relativistic effects of time delay, aberration and 
sweepback of the open field lines. While these effects are essential in under- 
standing pulse profiles, the phase-averaged flux is adequately described without 
a relativistic model. We use a polar cap acceleration model for the 7-ray emis- 
sion. We present the preliminary results of our recent study and the implications 
for observing millisecond pulsars with GLAST and AGILE. 



Subject headings: radiation mechanisms: non-thermal — magnetic fields — stars: 
neutron — pulsars: general — 7 rays: theory 
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Introduction 



Gamma-ray pulsars are the brightest objects in the sky above 100 MeV and the only 
identified GeV Galactic sources. Although there are presently around 1700 known radio 
pulsars, only a relatively small fraction of these are detected at higher energies. The EGRET 
detector on the Compton Gamma -Ray Observatory (CGRO) detected six 7-ray pulsars with 
high confidence (jThompsonll2004j ) . five of which were known radio pulsars. In addition, there 
were a few 7-ray pulsars detected with lower confidence, one of which was the millisecond 
PSR J0218+4232. The Large Area Space Telescope (GLAST), set for launch in late 2007, 
will have a sensitivity about 30 times better than that of EGRET and is expected to detect 
many more 7-ray pulsars. 

Population synthesis can predict the number of expected 7-ray pulsar detections, but the 
results are highly model-dependent. As a consequence of this sensitive model dependence, 
the number of 7-ray pulsars that GLAST does detect, especially the ratio of radio-loud to 
radio-quiet pulsars, will be an excellent discriminator between 7-ray mode ls. Using polar 



cap models to describe the 7-ray emission, population synthesis studies of iGonthier et al. 



( I2OO2I . I2OO4I ) have estimated that GLAS T should detect several hundred pulsars. Population 
synth esis studies for outer gap models (I Jiang fc ZhangI l2006l : iHarding. Grenier fc Gonthiei 
20061 ) generally predict fewer radio-loud 7-ray pulsars and a much higher ratio of radio-quiet 
to radio-loud pulsars. These studies evolved neutron stars in the Galaxy from birth distri- 
butions of period, magnetic field, position and space velocity, but excluded the population 
of millisecond pulsars (MSPs). Since MSPs are thought to be recycled through spin-up by 
accretion from a binary companion, their evolution is somewhat more complicated. 

Because of their very short periods, MSPs can have spin-down luminosities that are com- 
parable to those of young pulsars and therefore should be sources of high-energy photons. 
Indeed, a much higher fraction of the total radio MSP population (around 225) have been de- 
tected as X-ray sources (around 60, including Globular Cluster pulsars) compared to normal 
pulsars. The puls ed X-ray spectra of most MSPs in the Galactic plane are dominated by non- 
thermal emission (IKuiper et al.ll2003l ). so these sources must be accelerating particles to high 
energy. Polar cap models predict that MSPs should produce a high-energy emission compo- 
nent due to curvature radiation up t o energies around 50 GeV (jllarding. Usov. &: Musliniov 
2OO5I : ILou. Shibata. &: Melrosell2000l ). Outer gap models also predict that MSPs could pro- 
duce high-ener gy emission from part icles moving downward toward the stellar surface from 
the outer gap (jZhang fc Chengll2003r). Predictions for the nurn ber of 7-ray MSPs in globular 
clusters have been discussed by IWang. Jiang, fc Chengj (120051 ) . 



Studies of the radio characteristics of MSPs indicate some differences with those of 



the normal pulsar population. iKramer et al.l (119981 ) and also iBailes et al.l (119971 ) found that 
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MSPs are somewhat less luminous than normal pulsars. However, these studies compare the 
pseudo-luminosities, i.e. the observed radio flux times the square of the distance, which are 
dependent on solid angle and observer viewing angle, and do not address real luminosities. 
These studies also found that the beam radii of MSPs, although larger on average than 
those of normal pulsars, are smaller than expected by an extrapolation of beam radii with a 
P~^l'^ dependence of the normal pulsars population. The average MSP beam widths seem to 
reach a saturated value of around 30° — 50°, where as the extrapolated value is greater than 
100°. We will argue however that relativistic effects (IHarding. Grenier fc Gonthierll2006l ) may 
strongly influence the shape of MSP profiles. MSP pulse profiles have considerably larger 
duty cycles than no rmal pulsars and the beaming fractions are thought to be around 75% 
( iKramer et al.lll998l ). The profiles were originally thought to have more complex structure, a 
feature that was attribut ed to the presence of multipole field com ponents, possibly induced 
by the accret i on ph ase (IKrolik et al.l Il99ll : iKramer et al.l Il997l ) . But the observations of 
Kramer et al.l (Il998l ) do not support this claim and show that MSP profiles display in general 
only marginally more complexity. The more observable details afforded by the larger duty 
cycles make the profiles seem more complex with the large beams increasing the chance 
of detecting emission from both poles. Additionally, the emission height dependence on 
frequency (radius-to-frequency mapping) is less strong, and this is understandable since the 
MSP magnetospheres are much smaller. In general, MSP radio profiles are consistent with 
an expanded version of normal pulsar emission on open dipole field lines, with the smaller 
widths being attributed to emission occurring only on a subset of the open field lines or to 
relativistic distortions of the beam. 

We present in this paper new results of a population synthesis of MSPs born within 
the Galactic disk. In this study we define MSPs as those pulsars with period derivatives 
log(P) < —19.5 — 2.51og(P). This criterion includes the binary pulsars J1744-3922 and 
B0655+64 with periods of 172.4 ms and 195.7 ms, respectively, but excludes the binary 
pulsar J1711-4322 with a period of 102.6 ms. By means of a Monte Carlo population code, 
we treat them as point particles and evolve their trajectories, periods and period derivatives 
from their birth (at time of their last spin-up phase) forward in time to the present. We 
exclude MSPs in globular clusters from our study. After determining a present-day spatial 
equilibrium distribution, we assign radio and 7-ray characteristics to each MSP and then filter 
its properties for detection through a select group of radio surveys and the 7-ray instruments 
EGRET, AGILE and GLAST. We normalize the number of simulated radio MSPs to the 
number detected by the group of radio surveys, allowing us to predict the number of radio- 
loud and radio-quiet 7-ray MSPs detected by EGRET, AGILE and GLAST. 
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2. Selected Radio Surveys 

We recently added the Swinburne Intermediate Latitude survey to our selected group 
of radio surveys used in previous studies. Of the ten surveys, six of them, Arecibo 2 & 3, 
Greenbank 2 & 3, Molongo 2 and Parkes 2, have an observing frequency around 400 MHz, 
and four surveys at 1400 MHz include Parkes 1, Jodrell Bank 2, Parkes Multibeam, and 
Swinburne Intermediate Latitude. In our computer code, we use the characteristics of these 
ten surveys to determine the m inimum flux threshold for each simulated pulsar (for details 



see 



Gonthier et al.l (120021 . |2004J ) ) at the observing frequency of the survey. While we recognize 
that some of these radio surveys were not very sensitive to the detection of MSPs, we include 
them to confirm null detections. When the simulated radio flux is below the flux threshold 
Smin of all the radio surveys that potentially are viewing the simulated pulsar, we refer to 
the pulsar as being radio-quiet. Otherwise, the pulsar is radio-loud. These characterizations 
are relative to the sensitivity of the surveys and are not addressing the nature of the emission 
process. We then are restricted to comparing the characteristics of the pulsars simulated by 
our computer code to those MSPs detected by this select group of radio surveys, excluding 
MSPs detected only by other surveys. 



3. Kick Velocity and Equilibrium Spatial Distribution 



We a s sume exponential spatial distributions of the form described in the study of 
Paczyhskil (Il990l ). but using a scale height parameter above the Galactic plane of 200 pc 
for MSPs in contrast to 75 pc used in our simulations for NPs. We use a Maxwellian distri- 
bution with a characterist ic width of a ,, = 70 km/s to describe the supernova kick velocity 



at birth from the study of iHobbs et al.l (120051 ). resulting in an average velocity of 110 km/s. 
We assume a uniform birth rate back in time to 12 Gyr. To determine an equilibrium spatial 
distribution, we evolved a group of 4 0,000 neutron s t ars to the present without ma king any 



selections. As in our recent studies, iGonthier et al.l (|2007 ) and 



the Galactic potential (mass model 2) of iDehnen fc Binneyl (Il998l ) (W. Dehnen private com- 



Story et al.l (120061 ) . we use 



municat ion) to evolve the trajectories of MSPs using a 5th order Cash-Karp Runge-Kutta 
routine (jPress et al.lll992l ). While MSPs are born in binary systems, we have not attempted 
to simulate the evolution of the binary system, but rather treat the system as a point par- 
ticle and evolve it within the Galactic potential. With low kick velocities, the MSPs remain 
bound to the Galaxy and achieve normalized equilibrium distributions in Galactic radius 
(R) and Galactic height (Z) shown in Figure 1 as solid curves with the initial distributions 
represented by the dotted curves. 



Fitting the equilibrium distributions results in radial and out-of-plane scale heights of 4.2 
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kpc and 0.50 kp c, respectively. The Z scal e height of 0.50 kpc is in good agreement with the 
one obtained by lCordes &: Chernofa (119971 ) of 0.50 kpc ar id with the recent determinat i ons o f 
the scale heights of low mass X-ray binaries (LMXB) by I Grimm. Gifanov fc SunyaevI ((2002) 
of 0.41 kpc. In subsequent calculations, we use these equilibrium distributions to randomly 
select the present-day spatial distribution of MSPs within the Galaxy. In choosing the R and 
Z values of MSPs from such an evolved distribution rather than evolving each MSP from its 
birth location, we ignore any dependence of evolved R and Z on MSP age. 



According to IShklovskiil (1l970l ). the transverse motion of a pulsar causes a change in 
the pulsar's detected period. Due to the Doppler effect, an increase in the pul sar's period 
deriv ative is observed. The contribution to P due to this effect is given by ([Manchester 
1999h 

^^ y f 

pc 



p 



2.43 X 10"^^ ( — I , 

\ms / ymas/yr 



(1) 



where P is the pulsar period, fi is the proper motion, and d is the distance. For MSPs, whose 
period derivatives are extremely small, this effe ct can increase period derivatives significantly, 
in some cases by as much as 90% as noted by lToscano et al.l (119991 ) who have also included 
the smaller effects of Galactic rotation and vertical acceleration. In this study, we only 
include the Shklovskii effect due to the proper motion of the simulated pulsar. 



4. Magnetic Field and Initial Period 



While it may be that the magnetic field of MSPs decays slowly, we have assumed a 
constant field during their lifetime. However, due to selection effects, the observed field 
distribution will differ from the initial distribution. We explored power laws with various 
indice s to characterize the distr ibution of pulsar surface magnetic field at the birth line. 
While ICordes fc Chernofj (119971 ) preferred a power law with an index of -2, the group of 22 
MSPs used in their study had smaller magnetic fields. For the group of 56 MSPs in this 
study, we find a preferred index of -1 with a normalized distribution of the form 



n{Bs) 



1 



B. In 



Brr 



(2) 



where 5s is in units of 10*^ G. We choose a Br, 



1 and a Br, 



10^ that provides the 



best agreement between various simulated and detected distributions. 

Having the n aagnetic field for a simulated MS P, we can obtain its initial period with a 
spin-up relation. lAlpar et al. Jl982h and similarly iBhattacharya &: van den Heuvell (Il99ll ) 
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found that for recycled pulsars spun up by accretion at the Eddington limit from a binary 
companion, the initial period is estimated by the expression 



P. 



O.SPg^Vs, 



(3) 



which can be rewritten as a birth line in the P — P diagram by the form 

logPo = ^logPo-15.4, 



(4) 



where the period is in seconds and the period derivative is s ■ s^ . Having the magnetic field, 
the initial period, and the age, we can determine the present period and period derivative of 
the MSP assuming a pure dipole spin down. However, recent studies of LMXBs where the 
spin periods of accretin g neutron stars we re measured have allowed the estimation of the 
surface magnetic fields. iLamb fc Yul (120051 ) conclude that the general properties of LMXBs 
can be understood if the accretion rates range from the Eddington critical accretion rate 
to 5 X 10~^ times that rate. These different accretion rates lead to different birth lines in 
the P — P diagram. To incorporate a distribution of accretion rates, we have included an 
approximate procedure in which we dither the intercept of the following birth line 



log ^0=3 Po 



14.9 - (5, 



(5) 



where the dithering parameter 5 varies from to 2.8. The birth line can be reformulated by 
the expression 

Po = 0.18 X lO^'^/^Pg^Vs. (6) 



The extremes of the dithering parameter 5 represent the a pproximate Edding ton accretion 
critical rate and 5 x 10~^ times that rate (see figure 4 in iLamb fc Yul (120051 )). While we 
explored a Gaussian distribution of the dithering parameter, we obtain better agreement 
with a ramp distribution that increases by a factor of 4 between and 2.8. With the 
magnetic field randomly selected from the distribution of Equation 2, we then use Equation 
6 to obtain the initial period of the simulated MSP. However, we impo se a minimum initia l 
period of Po^^^ = 1.3 ms in accord with the RXTE studies of LMXBs by IChakrabartyl (120051 ). 
Having randomly determined the pulsar's age assuming a uniform birth rate, we spin the 
pulsar down to obtain the present period and period derivative assuming constant magnetic 
surface field. A constant magnetic field used in this study of MSPs is in contrast to the 



assum ption of magnetic field decay in our simulations of NPs (IGonthier et al.l 120021 . 12004 
20o3). 
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5. Radio Luminosity and Beam Geometry 



In our previous studies of NPs (iGonthier et al.ll2002l . 12004 120071 ). we adopted the in- 
trinsic radio luminosity model given by the general form 



(7) 



from lArzoumanian. Chernoff fc CordesI (120021 ) (ACC) with Lo = 2.lx lO^^m Jy • kpc^ • MHz 
a = —1.3 and = 0.4, but reduced the luminosity Lo by a factor of 60 fiGonthier et al. 



20041 ) and by 73 (IGonthier et al.ll2007l ) to obtain adequate agreement between the simulated 
birth rate, flux, and distance distributi ons and those detecte d, particularly by the Parkes 
Multibeam Pulsar survey (PMBPS) (see IGonthier et all (|2004j )). In ACC, the predicted birth 
rate was 0.13 neutron stars per century. In order to obtain a birth rate near 2 neutron stars 
per century, we had to decrease the luminosity constant Lo- The assumed radio luminosity 
is directly related to the simulated neutron star birth rate. Various studies propose neutron 
star birth rates that span somewhat of a range. For example, the a ran ge of birth rates 



from .9 to 1.9 per century was obtaine d from the analysis of PMBPS by IVranesevic et al. 



(I2OO4J ). In a recent population synthesis, iFaucher-Giguere &: Kaspil (120061 ) estimate a larger 



birth rate of about 2.8 per centur y, whereas 



Lorimer et al. 



J2006h ind 1.4 ± 0.2 per century 



from their analysis. More recently IGonthier et al.l (120071 ) determined the reduction factor of 
the radio luminosity by normalizing the simulated ne utron star birth rate to rate of Typ e 
II supernovae of 2.1 per century, following the work of iTammann. Loffer fc Schroderl (Il994j ). 
using only the Parkes Multibeam survey, as we may have the best description of the flux 
threshold Smin for that survey (Crawford, private communication). 

While the neutron star birth rate of 2.1 per century provides a constraint on the lumi- 
nosity of NPs, MSPs appear in a very different region in the P — P diagram, requiring careful 
consideration of the period and period derivative dependence of the radio luminosity used 
by the population synthesis study. Often, radio astronomers refer to the radio luminosity in 
terms of a "pseudoluminosity" {S(f) determined from the radio flux 5* observed at a specific 
frequency and the pulsar distance d. One cannot generally determine the intrinsic radio 
luminosity of a given pulsar from the Sd'^, even at a given frequency, as the detected average 
flux is strongly dependent on the viewing geometry. However, if the number of detected 
pulsars in a group is large, one might be justified in assuming that the emission region of 
the group of pulsars is completely sampled in a random fashion. Therefore, the detected 
Sd"^ would be proportional to the intrinsic luminosity of the group of pulsars if the surveys 
represented an unbiased sample of the true flux distribution. However, since radio surveys 
necessarily sample the high end of the flux distribution, inferring the intrinsic luminosities of 
pulsars is difficult at best. In addition, selection effects of the radio surveys may be different 
for NPs and MSPs, and these effects do influence the sampling of the emission region. For 
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example, the width of the pulse profile is used in the determination of the flux threshold 
Sinin for each radio survey; therefore, one would expect that narrow pulse profiles or larger 
inclination angles are preferably detected. Regardless of these ef fects, such compariso ns are 



important in characterizing differences between NPs and MSPs. Kramer et al.l (Il998l ) com- 
pared the average log{S(P) (mJy • kpc^) of a group of MSPs with that of a group of NPs 
and found that the mean \og{Sd?) values for their groups of 31 MSPs and of 369 NPs of 
0.5 ± 0.2 and 1.50 ± 0.04, respectively. The authors conclude that MSPs and NPs have 
very different ra dio luminosity distrib utions. Since this study was performed, the NE2001 



distance model (ICordes &: Lazidl2002l ) was developed. We use the pulsar distance given by 



the NE2001 model unless the distance has been determined by other more accurate methods. 



have a recorded uncertainty (. 


jorimer et al. 


1995: 


Kramer et al. 


1998: 


Lorimer et al. 


2000; 


Manchester et al. 


2001 


Hobbs et al. 


2004 




Kramer et al. 


2003: Faulkner et al. 


2005. 


20041; 


Morris et al. 


2002 


Stairs et al. 


2005h. we find average losiiSd?) values of 0.38 ± 0.15 for our 



select group of 39 MSPs (24 with S1400 fluxes) and 1.17±0.02 for a group of 1102 NPs (869 
with S1400 fluxes (references in ATNF catalogue)) detected by Parkes Multibeam and the 
Swinburne Intermediate Latitude surveys. The uncertainties are just the statistical error in 
the mean and do not take into account the error in the fluxes and distances. Subtle selection 
effects may be introduced within the context of a population statistics study with a group 
of surveys that have different sensitivities. As a result, it rnay be inappropriate to compare 
these mean \og{S(P) values with those from iKramer et al.l (119981 ) because of the different 
samples of pulsars. The MSPs recently discovered by the Parkes Multibeam and Swinburne 
surveys are not only more distant but also weaker with effectively smaller \og{Sd?'). We have 
chosen to reproduce the \og{Sd?) values of the normal pulsars and MSPs detected by these 
surveys to constrain our luminosity model. 

We adjust our intrinsic radio luminosity model such that the simulated mean \og{Sd?) 
values approximately match both the mean values of 0.38 ± 0.15 for MSPs and 1.17 ± 0.02 
for NPs detected by both the Parkes Multibeam and the Swinburne Intermediate Latitude 
surveys. We find that the exponents of the period (a) and period derivative (/3) dependence 
of the luminosity are not uniquely constrained to specific values but are correlated with each 
other. We find better agreement between the various simulated distributions for NPs with 
values of a = —1.05 and (3 = 0.37. The coefficient of the luminosity is adjusted to give a 
birth rate of 2.1 neutron stars per century using the Parkes Multibeam and the Swinburne 
Intermediate Latitude surveys. 

We were able to find a common radio luminosity that provides adequate results for both 
NPs and MSPs, given by the expression 

L = 1.76 X 10^°p-^ °^P°-^Wy ■ kpc^ • MHz. (8) 



- 9 - 



In this study, we assume that the radio emission originates from a single core beam 
and a single cone beam axially aligned with the magn etic field from each p ole of the star 
and isotropically distributed in the sky. As indicated by iKramer et al.l (119981 ) , the profiles of 
MSPs display similar structure to those of NPs, suggesting that the core or central beam and 
the conal beam make similar contributions to NP and MSP profiles. However, as discussed 
later, due to the rapid rotation of M SPs, relativistic effe cts are more evident in the pulse 
profiles of MSPs than in those of NPs (IXilouris et al.lll998r). W e use the same description for 
the core beam as in our pa st study of Gonthier et al.l ( 2004 ). but we replace t he geometry 
of the cone beam from the Mitra &: Deshpandel Jl999h with the description of iKijak &: Gill 
(Il998l . booah having the form 



Pc 



(9) 



which describes the characteristic width of the conal beam at 0.1% of the peak intensity of 
the profile assuming that the edges of the pulse are along the last open field line. They find 
that the altitude of the emission tkg in stellar radii is given by 



(10) 



where P-15 is in units of 10 s • s . In this formulation, the overall period dependence 
of the characteristic width of the conal beam is / ?r,nne ^ -P "'^^/ reproducing well the trend 



observed for NPs and MSPs (IKramer et al.lll998l : lKramerll2002l ) 



In the formulation of a Gaussian conal beam, the characteristic width Pcone requires the 
specifica tion of the radius ^ and width We of the annulus of the beam. The characteristic 
width of lKijak fc Gill (Il998l . l2003l ) is for 0.1% of the peak intensity making pcone = ^+2.63we. 
Within such a model, one is free to specify some form for either 9 or We- We adopted 
We = O.lSpcone in our study of both normal and MSPs. We find that the population statistics 
are not very sensitive to the choice of this value, which we have arrived at from our study of 
three peak pulse profiles of radio pulsars. In our study of the characteristics of core and cone 



beams of about 20 radio pulsars whose profiles have three peaks, we found (iGonthier et al. 



(I2OO6I ) and in preparation) that the ratio of core-to-cone peak fiuxes does not follow the 
dependence suggested in the work of ACC. In the ACC model, the profiles of MSPs 
are dominated by the central core beam, so tha t even at sma l l imp act angles the conal 
beam would not be seen. As seen in the study of IKramer et al.l (119981 ). similar complexity 
is manifested in the MSPs profiles as in profiles NPs, indicating that the conal beam is 
quite prominent. Recent pol arimetry studies of young radio pulsars also call the ACC core- 
to- con e ratio into question. ICrawford. Manchester fc Kaspil (I2OOII ) and I Crawford fc Keim 
(I2OO3I ) measured the radio polarization of nine pulsars, six of which had characteristic ages 
less than 100 kyr. They found that the profiles of all the young pulsars indicate a significant 
degree of linear polarization and a low degree of circular polarization, which is a traditional 
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feature of conal emission. In this study, the two older pulsars with characteristic ages greater 
than 1 Myr do not show any degree of polarization. They conclude, as did [Manchester 



( Il996l ). that the profiles of young short period pulsars tend to be c haracterized by partial 
cones. Similarly, in a more recent study [Johnston fc WeisbergI (120061 ) measured polarization 
profiles of 14 young pulsars with ages < 75 kyr. They found that the profiles are dominated 
by linear polarization, suggesting that the core beam is weakly manifested, if at all. Because 
of these inconsis t encies , we have reassessed the relationship between the core and cone beams 
( iGonthier et al.l (120061 )). using three-peak normal pulsars and a few three-peak MSPs. In 
our three-peak pulsar study, we find for the ratio of peak fluxes a power-law dependence 
with a break in period at about 0.7 s having the approximate form 



25Pi-3 

4P-1.8 



^Gr,^<0.7s 
z.-°f,P>0.7s. 



(11) 



Hence in this model, short period pulsars are much less core dominated than in the 
ACC model. The model do es not distingu i sh bet ween young short-period pulsars and MSPs, 
although there is evidence (iKramer et al. (Il999h) that the ratio of core-to-cone peak fluxes 
may be different for a young pulsar and a MSP with the same period. The shapes of 
MSP profiles are complicated by relativistic effects, but due to the wide emission beams, 
contributions from both poles are often seen in the profiles, making the analysis of the 
viewing geometry and conal widths difficult. With this model, we do accommodate more 
complex pulse profiles for short period pulsars than in the ACC model. 

In describing the core and cone angle-integrated spectra, we assume indices of -2.36 and 
-1.72, respectively. While these indices are not measurable, they are related in a complex 
fashion, due to selection effects, to the measured spectral indices. Using these values, our 
simulation gives a flux-averaged spectral index of —1.80 ± 0.05 for MSPs. Within our select 
group of surveys, we fl nd 20 MSPs within the ATNF catalogue whose spectral indices have 
recorded unce r tainti es (iToscano et al.lll998l : lLorimer et al.lll995l ) with a mean of —1.89 ±0.10. 
Kramer et al.l (119981 ) found a mean spectral index of — 1.8±0.1 for a group of 32 MSPs. They 
also point out that MSPs have signiflcantly steeper mean spectral indices than NPs. For 
a group of 346 NPs, they flnd a mean spectral index of —1.60 ± 0.04. In our simulation, 
we found it necessary to have softer angle-integrated spectral indices for the core and cone 
beams of NPs, resulting in a mean spectral index of —1.58 ± 0.01. This value agrees well 
with the mean spectral index of — 1.66±0.04 for the 245 NPs with indices that have recorded 
uncertainties in the ATNF catalogue (references therein). 

Adopting this radio luminosity and beam geometry model, we attempt to describe radio 
emiss ion properties from both NPs and MSPs (IGonthier et al.ll2007l : iHarding. Grenier fc Gonthier 
20061 ) . The success of this radio model is remarkable in describing the physical luminosities of 
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both NPs and MSPs, given the range of three decades in period and seven decades in period 
derivative between NPs and MSPs. However, it is clear that there are many complexities 
that may not modeled correctly and that the population synthesis study may not be very 
sensitive to various characteristics of the beam geometry. 



Relativistic Effects in the Radio Profiles 



Because of their rapid rotation, MSPs can display very different features in their pulse 
profiles than those observed in the profiles of NPs. These features can be partly understood 
in terms of t he relativistic effects of aberration, time delay and the s weep back of the retarded 
dipole field (jPyks fc Harding] |2004J : iDyks. Harding fc Rudakll2004l ). The shorter the pulsar 
period, the more dramatically these effects are manifested in the profiles. According to 
Equation 10, the conal emission of pulsars with short periods occurs at high altitude relative 
to the light cylinder radius. For MSPs the altitude of emission for the cone beams {tkg) is 
~ 0.2 — 0.6 Rlc- The combined effects of time delay and aberration result in a shift in phase 
of the conal beam. Deformations in the shape of some prof iles of MSPs are a resu lt of the 
distortion of the open field volume of the retarded dipole (jPyks &: Harding |2004| ) . While 
relativistic effects are important in understanding the s hape of the profil e and in studying 
correlation with high-energy profiles, in a recent study (IStory et al.ll2006l ). we find that the 
phase-averaged flux is very similar for relativistic and nonrelativistic emission. Since in 
population statistics studies, the average flux of the profile is the determining parameter 
that is compared to the threshold of each radio survey, we have neglected these effects in 
this study. 



7. Gamma-ray Luminosity and Beam Geometry 

Because of the lack of a theoretical understanding of the mechanism that produces the 
radio emission, the radio emission model is phenomenological by nature. On the other hand, 
the 7-ray emission models are motivated by theory. Two competing 7-ray emission models, 
the polar cap and outer gap models, assume different locations in the magnetosphere where 
the acceleration of charge occurs. In the outer gap model, the acceleration of charge takes 
place in vacuum gaps formed along the last open field line, between the null surface and 
the light cylinder. While we are be ginning to incorporate an outer gap in our simulations 



(IHarding. Grenier fc Gonthierll2006l ). treating it on the same footing as the polar cap model 



using the same group of evolved neutron stars and selection criteria, we use only the polar 



cap model in this study of MSPs. Recent outer gap models for MSPs (IZhang &: Chengll2003l ) 
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actually assume that the high-energy emission originates not from the outer gap but near 
the neutron star, from particles flowing down from the gaps, and thus would have a very 
different high-energy emission geometry from that of normal pulsars. 

Two regions in the P — P diagram, separated by the curvature radiation (CR) pair death 
line, differentiate the 7-ray luminosities and beam geometries within the polar cap model. 
Pulsars above the CR line can produce CR photons that are energetic enough to pair produce 
in the strong magnetic field above the neutron star surface. The CR-initiated pair cascades 
have sufficient multiplicity to fully screen the parallel electric field everywhere except in 
the slot gap along a narrow region near the last open field line where the primary particles 
are a c celera ted along unscreened electric field lines. Recent studies ( iMuslimov &: Harding 
2OO3I . I2OO4I ) of the emission above the CR pair death line in the polar cap model have 
provided a framework for understanding 7-ray emission from the slot gap. The electrons 
in the low alti tude slot gap initiate pair c ascade emissio n, which was incorporated i n our 



simulations of iGonthier et al.l (12004 . 



20071 ). In addition, iMushmov fc Hardind (120041 1 find 



that the primary particles along the last open field lines can be accelerated to high altitudes as 
the potential in the s lot gap remains unscre ened and the emission beam forms a characteristic 
caustic component (jPyks fc RudakI 120031 ). In recent population stu dies of NPs, we have 
begun to include both the low and high altitude emission geometries (IGonthier et al.l 120071 : 
Harding. Grenier &: Gonthierll2006l ) and both are included in this study for those simulated 
MSPs above the CR line. 

However, as one can see in Figure 2, all the MSPs with the exception of one lie below 
the CR line where curvature radiation no longer produces pairs. These pulsars can produce 
much weaker pair cascades through inverse Compton scattering radiation. When pulsars 
are no longer able to produce abundant pair cascades, they do not form slot gaps above 
the polar cap iMuslimov fc Harding] (120031 ) as the parallel electric field is not screened on all 
open magnetic field lines. As the 7-ray energies are below the pair threshold, the curvature 
radiation from the primary electrons escapes out to infinity along all open field lines. In 
this pair-starved polar cap model, the primary particles continue to accelerate and radiate 
to high altitudes above the polar cap out to the light cylinder along all the open field lines. 
However, CR emit ted by the particles can be inhib ited by resonant cyclotron absorption 
of radio emission (IHarding. Usov. fc Muslimovl |2005| ) where synchrotron radiation is more 
efficient particle energy loss mechanism than CR. Since synchrotron photons usually have 
lower energies than curvature photons, we have not included this mechanism in the pair- 
starved model in this study. The CR em ission model below the CR p a ir dea th line used in 
this study is developed from the work by lHarding. Usov. fc Muslimovl (120051 ). 



The inclination angle a and the viewing angle C, determine which open field lines are 
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sampled by the line of sight. As illustrated in Figure 3, the line of sight intersects tangentially 
a particular open field line with a radius of curvature pc at a radial distance r and polar angle 
9. The open field lines are defined within the polar cap angle given by the expression 



pc 



sm 



2^R 



1/2 



(12) 



where P is the pulsar period in seconds, R is the neutron star radius (10^ cm) and c is the 
speed of light. We partition the open field lines through a dimensionless parameter, ^, that 
varies between and 1 and is defined as 



9. 



(13) 



pc 



where 6s is the polar angle of the open field line at the intersection with the surface of the 
star. The particle emits curvature photons along the line of sight at an angle 9^ = 39/2 
given by 

cos 9^ = sin a sin ( cos (j) + cos a cos (, (14) 
where is the phase angle, a is the magnetic inclination and ( is the viewing angle. 



W e approximate the accelerating electric field from Equation 14 of iHarding fc Muslimov 
( 1l998l ) using the expression 



sm 
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;i-^^) statvolt-cm~\ (15) 



where k = 0.15 is the general relativistic inertial frame dragging factor, B12 is the surface 
magnetic field in units of 10^^ G, rj =r/R is the dimensionless radius and if is the azimuth 
angle around the magnetic pole given by the expression 

cos ( — COS a cos 9^ 



cosy) 



sin a sin 9-, 



(16) 



The gain in energy of the accelerating primary electron is compensated by the CR losses 
faarding Muslimov fc ZhaneJbooi ISuhk. Rudak fc Dvksl[200ol : ILou. Shibata fc Melroselboooh , 
so that the electron Lorentz factor 7 becomes radiation-reaction limited to 



3p^ 



e Ell 



(17) 



where pc is the classical radius of curvature of the field line given by 

r (1 + cos^ 9f^^ 



pc 



3sin^(l + cos^)' 



-14- 



At low altitudes near the stellar surf ace, r < R(sm9pr./3 + 1), th e near surface electric field 
is approximated from equation 21 in lHarding &: Muslimovl (120011 ) by the expression 



E|| = -4.05Si2sin3^p,(r7-l) 



0.15(1 -e 



2.19N0.705 



cosa + -sin^p,e " (1 



-2N0.65 



sm a cos cp 



(19) 

and goes to zero at the stellar surface. The CR spectrum therefore has a spectral index of 
-2/3 with a high-energy cutoff ecR given by 



(20) 



in me(? units where e is the electron charge and Ac = h/ {mec) is the electron Compton 
wavelength. 



The instantaneous curvature radiation power spectrum has the form (jJacksonlll975l ) 



^CR 



(21) 



where e is the photon energy, 7 is the particle Lorentz factor, ctfine is the fine structure 
constant, and the k(x) function is defined as 



K[X) = X 



K,/,{x')dx' 



22/3r(|) X^^X^l, 

1.253x^/V^,x > 1 



(22) 



We find that the low energy asymptotic form provides a sufficiently accurate description of 
the power spectrum above e^=100 MeV and the instantaneous curvature photon radiation 
spectrum per primary is given by 
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(23) 



which is similar to equation 7 in iHarding. Usov. fc Muslimovl (120051 ) . This expression can 
be integrated from e^=100 MeV to the high-energy cut off ecR to give the expression 



NcRi> e^,e^,^,Lp) 
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(24) 



which provides the number of curvature photons emitted along the line of sight 6.y from an 
open field line ^ at a radial distance r per primary particle, where r is given by the expression 



sin^(^6'pc) 



(25) 



- 15 - 



The Goldreich- Julian current of primary particles from the polar cap is uniformly distributed 
over the polar cap and is given by the expression 

iVp = 1.3 X 10^°5i2P"^ particles/s. (26) 

The number of primary particles in a particular patch on the stellar surface at R, 6's, and 
which gives the number of primaries along a particular field line ^, is provided by the 
expression 

^GJ = 7m FT- ^^^> 

Itx[ \ — cos Upc) 

The total photon luminosity from one pole can be obtained by integrating over the open 
field volume to give 

-2tv fl i-Rlc 



'>s^)=hGjl I I Ncr{> e^,e^,^,ip)epcsm^epcdrd^d(p. (28) 
Jo Jo Jr 

dr\ 2R sin 6 cos 9 , , 
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With 

{^dQ) sm.\iep, 
the integral over r can be written as an integral over Q to give 

T \ r f' r""""^ ( 2RhGjNcR{> e,,e„^,^)ep,cose \ . 

Aotai(> £j) = / / mTTl smOded^dcp. (30) 

Jo Jo Je{R) \ csm(^^pj / 

where the integral over ^ (the Opc has been included in the integrand) is over the open field 
lines and ^9pc is the angle at the footpoint of the field line on the stellar surface. The 
expression in parentheses represents the photon emission from a particular field line at ^ and 
?7, integrated above e^=100 MeV, along the line of sight into a solid angle dVL and is given 
by the expression 

dNcR{> e^, ^) _ ^■^5'^O^f^neNGJCR9p,p-^/^ COS 9 (e'f - 

dh ~ COS Op,) (hey sm{^9pc) ' ^^^^ 

The total emission from the polar cap for a particular phase bin is obtained by integrating 
over parameter ^ to include the contributions of all the open field lines along a line of sight 
defined by the viewing angle (. The total emission in the phase bin is given by 

^„(^^^_^)_^'*W>|^«^,j_ (32) 

Having calculated the 7-ray pulse profile for each MSP for a given viewing geometry, we 
average the profile to obtain the average photon flux, and compare it to the instrument 
threshold. Our resulting 7-ray efficiency (7%) for the total spectrum a nd profiles for the 
nearby MSP PSR J0437+4715 are in agreement with those obtained by IVenter fc DeJager 
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8. Gamma-ray All Sky Threshold Maps 

From the simulated 7-ray pulse profile, we obtain an average flux that we compare to all 
sky threshold maps for EGRET, AGILE and GLAST Large A rea Telescope (LAT). We us e 



the recently revised EGRET map that includes the dark clouds (ICasandjian fc Grenierll2007l ) . 
which has led to a radical reassessment of the EGRET unidentifled sources. The GLAST 
threshold has been improved and updated (Grenier & Casandjian private communication) 
as a 1 year GLAST LAT threshold map. The all sky map for AGILE (Pellizzoni private 
communication) has not been recently updated in our computer code. The detection of radio 
and 7-ray point sources within the code are independent of each other, allowing the tagging 
of radio-quiet (below the survey flux thresholds) and radio-loud 7-ray MSPs. 



9. Results 

To improve the simulated statistics, we run the simulation to obtain ten times the 
number of detected MSPs and then normalize the distributions accordingly. In Figure 4, we 
present Aitoff projections for the detected (a) and simulated (b) MSPs. Since MSPs are closer 
to us and are much older than NPs, the graphs indicate larger out-of-plane distributions than 
those of NPs. 

Of the ten radio surveys that we have included in our simulation, the surveys most 
sensitive to MSP detection are Parkes 2, Parkes Multibeam and Swinburne Intermediate 
Latitude which together account for majority of the 56 MSPs, as seen in Table 1. Within 
the limited statistics, there is good agreement between the number of MSPs detected and 
simulated among the ten surveys used in our selected group. The few detections in radio 
surveys that are not very sensitive to the detection of MSPs are also well reproduced by the 
simulation. 

The distributions in the P — P diagram are compared in Figure 5 for detected (a) and 
simulated (b) MSPs. The dotted broken lines represent the pair death lines for curvature 
radiation (CR) and for nonresonant inverse Compton scattering (NRICS). Below the CR 
death line, CR no longer is able to produce pairs. However, a limited number of pairs can 
still be produced via inverse Compton scattering above the NRICS death line. Below the 
NRICS death line, pair production is no longer possible, presumably inhibiting radio emission 
mechanisms. The upper and lower MSP birth lines indicated by the solid lines represent the 
approximate Eddington critical accretion rate (upper line) and 5 x 10~^ times that rate 
(lower line) discussed earlier. With the magnetic fleld distribution and the uniform 

distribution of birth lines between the upper and lower MSP birth lines, we reproduce the 
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observed broadness of the pulsar distribution in the P — P diagram. 

In Figure 6, we present histograms of the indicated measured and derived pulsar char- 
acteristics for the detected (shaded histograms), and the simulated (unshaded histograms) 
MSPs. Although the poor statistics make it difficult to judge the agreement, the simula- 
tion appears to reproduce the observed distributions fairly well. The clump of detected 
MSPs observed in Figure 5a is exhibited as narrow peaks in the period and period deriva- 
tive distributions. The narrow peak in the period derivative distribution is not present in 
the simulated distribution. In the comparisons of the age distributions, we calculated the 
traditional characteristic age P/ (2P) for the simulated pulsars. In the comparisons of the 
distances and d i spersi on measures of the MSPs, we have used the new distance model of 



Cordes &: Lazid (120021 ). The distance of the detected pulsar is established typically by its 
measured dispersion measure and location in the sky. The new distance model is used to 
determine the simulated dispersion measure from the simulated distance and location in the 
sky for each MSPs, which is then used in the calculation of the flux threshold Smin for each 
radio survey. Both the simulated dispersion measure and distance distributions indicate good 
agreement with the corresponding distributions of detected MSPs with a few more simulated 
closer pulsars. The simulated radial and height distributions within the Galaxy are consis- 
tent with those distributions of the detected pulsars. However, the simulated spectral index 
distribution is narrower than that of the detected pulsars. For the pulsar characteristics in 
Figure 6, we have performed a Kolmogorov-Smirnov (KS) statistics test of the histograms 
assuming a significance level (a) of 5%. We indicate in Figure 6 the resulting p- values for 
each of the comparisons. Sometimes the D statistic and the critical value are quoted in the 
literature and compared to each other to decide on the null hypothesis. The critical value is 
derived from a and the number of samples. Equivalent comparisons can be made between 
the p- value and a. If the p- value statistic is larger than a, the null hypothesis is not excluded 
at the level of significance determined by a. Our results indicate that in all comparisons, 
the distributions of simulated pulsars are consistent with those detected with the exception 
of the distributions of the Galactic height z and S400 whose p-value is less that 0.005. In 
the case of the S400 distributions, the statistics are insufficient in the number of measured 
fiux densities at 400 MHz for an adequate comparison. We simulate slightly more pulsars 
with larger S1400 fiuxes. We applied the two dimensional KS test to t he un binned period 



and period derivative distributions using the prescription of iPress et al.l (120021 ) and obtained 
a p-value of only of 0.03 at the same significance level of 5% indicating that the detected 
and simulated P — P distributions are not entirely consistent with each other. While the 
KS tests might be helpful in such comparisons, the numerous assumptions within the sim- 
ulation, such as the radio luminosity, the pulsar distance, and birth distributions, contain 
systematic uncertainties. Therefore, we do not want to overemphasize the importance of 
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these statistical tests. 

The detection of a radio pulsar in the computer code is independent of the detection 
of a 7-ray pulsar, allowing the identification of radio-quiet and radio-loud 7-ray pulsars. By 
radio-quiet, we mean that the radio flux of the simulated pulsars is below the flux threshold 
of the surveys that could potentially observe it; otherwise, the pulsar is flagged as radio- 
loud. In Figure 7, we show the simulated radio-quiet and radio-loud 7-ray MSPs for the 
instruments EGRET (upper right), AGILE (lower left) and GLAST LAT (lower right). For 
comparison, we also show the 7-ray pulsars that were detected by EGRET (upper left), 
including the MSP J0218-I-4232. In the case of Geminga (cross in upper right), it is not 
clear whether its radio silence results from a misalignment of the 7-ray and radio beams or 
from its intrinsic radio weakness. While it is unclear how many of the EGRET unidentified 
sources might be MSPs detected as point sources, the simulation predicts that a few may be 
radio-quiet 7-ray MSPs. 

All the simulated 7-ray pulsars that are detected by EGRET, AGILE and GLAST lie 
below the curvature radiation pair death line in the P — P diagram (Figure 7). As discussed 
above, in this pair-starved pola r cap model, curvature radia t ion es capes along all open field 



lines. However, as indicated by lHarding. Usov. fc Muslimovl (120051 ) . curvature radiation can 



be limited when the accelerated particles undergo resonant cyclotron absorption of radio 
photons under special circumstances. The excited particles emit synchrotron radiation typi- 
cally at lower energies than curvature radiation. While we are in the process of developing a 
full three-dimensional study of these processes, in the present calculation we have not limited 
the curvature radiation in altitude. Therefore, our results may be somewhat optimistic. In 
Figure 7, the 7-ray MSPs cluster in the lower left of the P — P diagram because the 7-ray 
luminosity increases with decreasing P. 

In Table 2, we present the 7-ray MSPs statistics for radio- loud and radio-quiet pulsars. 
Due to the uncertainties in the radio luminosity, we vary the coefficient of the radio lumi- 
nosity in Equation 8 by ±10% to obtain a range in radio-quiet 7-ray pulsars. Since the 
simulation is normalized to the detected number of radio pulsars, there is no variation in 
the number of radio-loud 7-ray pulsars. Of the 33-40 radio-quiet 7-ray MSPs predicted to 
be detected by GLAST, about 6 have fluxes above the 10~^ photons/cm^/s level (Grenier, 
private communication), that GLAST may be able to detect through blind period searches. 
However, since 70% of the MSPs in the ATNF catalogue are in binaries, many of these radio- 
quiet 7-ray pulsars are expected to be in binary systems. It will be difficult for GLAST to 
detect pulsed emission from these MSPs in blind period searches due to timing uncertainties 
associated with the orbital motion. Without radio ephemerides, GLAST will detect binary 
MSPs only as point sources. 
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In order to explore the nature of the radio-quietness of the GLAST pulsars, we present 
in Figure 8 histograms of the radio flux at 1400 MHz of the detected (solid dark gray 
histograms) and simulated (open thin histograms) radio pulsars in (a) and the radio flux of 
GLAST radio-quiet (solid light gray histograms) and radio-loud (open thick histograms) in 
(b). As expected, GLAST radio-loud pulsars have large radio fluxes. However, the GLAST 
radio-quiet pulsars form two groups, one with a flux distribution very similar to the radio 
detected and simulated pulsars in Figure 8a and the other with a much lower flux distribution 
that would be difficult to detect without deep radio pointings. 

The Parkes Multibeam and Swinburne Intermediate Latitude surveys are the two most 
sensitive surveys in our select group of radio surveys. However, they cover a rather narrow 
strip along the Galactic disk, and MSPs seem to be rather scattered in space and are not 
as correlated with the Galactic disk as are NPs. As a result, many of the GLAST detected 
pulsars are radio-quiet because they happen to be outside the sensitive survey regions or 
have radio fluxes that are slightly lower than the survey thresholds. 

In Figure 9, we plot the impact angle f3 = ( — a a.s a function of the radio flux at 1400 
MHz for the simulated group of GLAST radio-quiet pulsars. We can see that the two groups 
of fluxes correspond to two groups of impact angles. The MSPs in the low flux group with 
large impact angles will not be observed as radio sources without deep radio observations. 
However, the simulation suggests that they can be seen as 7-ray point sources. In the pair- 
starved model of 7-ray emission, curvature radiation occurs along all open field lines all the 
way out to the light cylinder. While the intensity of curvature radiation decreases at higher 
altitudes for nearby pulsars, we predict that GLAST should be able to detect them as point 
sources, yet their radio emission will be very difficult to detect as the radio beam and the 7- 
ray beam are separated by large angles. However, the higher flux group with smaller impact 
angles would be detected by GLAST as point sources and potentially are detectable with 
pointed radio observations. 

As a result, our simulations suggest that it will be important to follow up GLAST 
detections with radio observations in order to discover a greater number of radio-loud 7-ray 
pulsars. The predicted energy spectra of MSPs are very different than those of either NPs or 
AGNs. In this pair-starved polar cap model, the energy cutoffs of the MSP 7-ray sp ectra are 



predicted to be in the 10-50 GeV energy region (IHarding. Usov. fc Muslimovll2005l ) suitable 
for GLAST. 

Clearly the number of detected radio-quiet and radio-loud 7-ray MSPs depends strongly 
on the detection thresholds of radio surveys and 7-ray instruments. In Figure 10, we present 
the Log(N)-Log(S) for the radio flux at 400 MHz (a) and 1400 MHz (b) and for EGRET (c) 
and GLAST (d). The light shaded histograms represent the simulated undetected source 



- 20 - 



counts. The dark histograms represent the simulated detected source counts by the select 
group of ten radio surveys at 400 MHz (a) and 1400 MHz (b). For EGRET (c) and GLAST 
(d), the dark and medium dark histograms represent the simulated source counts for radio- 
loud and radio-quiet 7-ray MSPs, respectively. 

The appearance of fractions of a pulsar in the figure is a result of the re-normalizing the 
distributions because the simulation was run for ten times the number of detected MSPs to 
improve the simulated statistics. Assuming a uniform distribution of sources in space with 
isotropic emission of equal flux, one expects to a slope in the Log(N)-Log(S) diagram of -3/2. 
Indeed the tails of these distributions have slopes very similar to -3/2. While for EGRET 
and GLAST the Log(N)s of simulated detections have a sharp decrease with decreasing flux, 
the Log(N)s of simulated radio detections have a more gentle roll over, which may reflect the 
larger number of factors that determine the minimum flux threshold for the radio surveys. 
Both the number of undetected radio and 7-ray sources follow a simple power law far below 
the sensitivity of the radio surveys and 7-ray instruments. 



10. Discussion 



We have presented a Monte-Carlo simulation of the Galactic-plane population of MSPs 
and corresponding predictions for the numbers of both radio-quiet and radio-loud MSPs that 
are detectable as 7-ray sources. In order to accomplish this, we made significant improve- 
ments in our population statisti cs Monte Carlo code by a dding the more reahstic descrip- 
tion of the Galactic potential of iDehnen fc Binneyl (119981) with a rnore a ccurate 5th order 
Cash-Karp Runge-Kutta trajectory i ntegration routine (IPress et al.lll992l ). using improved 
all-sky threshold maps for EGRET (jCasandjian fc Grenierl 120071 ) and GLAST (Grenier & 
Casandjian, private communicatio n) , and including the Swinburne Intermediate Latitude 
radio survey (lEdwards et al.ll200ll ). We have modified the intrinsic radio luminosity model 
of ACC to describe the total radio lum inosity of both normal and MSPs. From our study of 
radio pulsars with three-peak profiles (iGonthier et al.ll2006l . 120071 ). we find that the ratio of 
the core-to-cone peak fluxes for short period pulsars is much smaller for short period pulsars 
than the ratio predicted by ACC model. Although MSPs have on average larger pulse widths 
than NPs, this resulting radio beam geometry model provides a reasonable agreement to the 
average pseudoluminosities and pulse widths at 1400 MHz and pulse widths of both NPs and 
MSPs. We include a 7-ray beam geometry and luminosity mo del of the curvature radiation 
from all open field lines within a pair-starved polar cap model (jHarding. Usov. &: Muslimov 
20051 ) that accounts for 7-ray emission below the curvature radiation pair death line, which 
is applicable to MSPs. These improvements in our computer code have allowed us to study 
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the population statistics of radio and 7-ray MSPs to complement our previous studies of 
NPs from the Galactic disk. As a first order study, we believe that we have improved the 
simulation by including a beam geometry for MSPs and an intrinsic radio luminosity whose 
functional form allows for the description of the population statistics of normal pulsars. The 
actual beam geometry of MSPs may turn out to be signific antly differen t and require a dif- 
ferent model than the one proposed here from the work of iKijak fc Gill (119981 . |2003| ). The 
population statistics study is not sensitive enough to discriminate between beam geometry 
models as the average flux of the profile is the quantity that is compared to the S'^mS of 
the radio surveys. A detailed study including relativistic effects will be necessary along with 
high quality polarization data may enable establishment of systematic trends that allows for 
the development of an adequate beam model for MSPs. 

In this study, we do not attempt to describe the numerous MSPs within globular clus- 
ters. Limiting our study to the MSPs born in the Galactic disk, we treat them as point 
particles, evolving their location within the Galactic potential and their spin-down, after 
spin-up by mass accretion from a binary companion star has ended. This group of old, 
short period pulsars with low magnetic fields are given radio and 7-ray beam and luminosity 
characteristics and filtered through a set of ten radio surveys and the 7-ray instruments 
EGRET, AGILE and GLAST. Our simulation is run until the number of simulated radio 
pulsars is equal to the number detected by the same group of r adio surveys. The radio lumi- 
nosit y is adjusted in our study of NPs (IGonthier et al.l 120071 : iHarding. Grenier fc Gonthier 
20061 ) to give a birth rate of Type I I supe rnovae of 2.1 per century to agree with the 
study of iTammann. Loffer fc Schroderl (119941 ) and to reproduce the detected mean pseudo- 
luminosities and spectral indices of MSPs and NPs, after which there are no further adjust- 
ments performed in the simulation of MSPs. Allowing for a ± 10% variation in the radio lu- 
minosity, we predic t a birth rate of 4 — 5 x 10"*^ MSPs per centur y, which is in good agreement 
with t he studies of iLorimen (120051 ) (3 x 10"*^ per century) and iFerrario fc Wickramasinghe 
( I2OO6I ) (3.2 X 10~^ per century). While the population statistics studies of LMXBs contain 
numerous uncertainties, our simulated bi rth rate agrees with t he estimate of the birth rate of 
LMXBs from the preferred model (F) of Kiel fc HurlevI (l2006l ) (6.5 x 10~^ per century). We 
are also in agreement with lKramer et al.l (119981 ) as our description of the spectral properties 
of the core and cone beams for MSPs requires steeper angle-integrated spectral indices than 
those for the simulated NPs. 



In outer gap accelerator models (ICheng. Ho fc Rudernianlll986f) . MSPs are not expected 
to produce 7-ray emission above 1-2 GeV Jzhang fc Chenej j2o'o3i r Hirotani private commu- 
nication 2005) due to the smaller electric field and larger radius of curvature in the gap. 
Thus , detected emiss i on fro m MSPs above 5 GeV would imply an origin near the polar 
cap. IZhang &: Cheng (120031 ) propose that particles accelerated in the outer gaps could be 
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the origin of MSP radiation at super-GeV energies, but only if the cascade emission from 
the particles as they travel down toward the neutron star surface is visible. However, it 
seems unlikely that both the upward-going radio emission and downward-going 7-ray emis- 
sion could be visible to the same observer. In that event, one would expect the radio and 
high-energy peaks to be out of phase, and X-ray and radio peaks are often in phase for MSPs 
(e.g. B1821-24, J0437-4715, B1921+37). In any case, GLAST observations of MSPs will be 
an excellent model discriminator. 

Our simulations predict the number of radio-loud and radio-quiet 7-ray pulsars observed 
by the instruments EGRET, AGILE and GLAST. EGRET was only able to detect one 
MSP, PSR J0218-I-4232, during its nine-year mission. EGRET did not have the capability of 
performing blind period searches of pulsars and required reliable radio ephemerides to detect 
7-ray pulsations. However, these observations were prior to the Parkes Multibeam Pulsar 
Survey and the Swinburne Intermediate Latitude Survey that tripled the number of detected 
radio MSPs in Galaxy (not counting those in globular clusters). We predict that EGRET 
should have detected as point sources eight radio-loud and ~ 20 radio-quiet 7-ray MSPs, 
some of which could be associated with EGRET unidenti fied sources. However, the re cently 



revised all-sky EGRET threshold map with dark clouds (jCasandjian fc Grenien 120071 ). used 
in our simulations, has significantly affected the catalogue of EGRET unidentified sources. 
As the all-sky threshold map for AGILE has not been recently updated, our predictions for 
AGILE might be somewhat optimistic. However, the trend in the statistics for EGRET, 
AGILE and GLAST suggests that the number of radio-quiet and radio-loud 7-ray MSPs 
predicted to be detected by AGILE might be reasonable. The Shklovskii effect is a significant, 
positive contribution to the intinsic P of the pulsar resulting in much larger measured Ps 
and needs to be taken into account in population statistics studies. The effect of having 
to decrease the intrinsic P to accomodate this effect results is pulsars with lower 7-ray 
luminosities. We predict ~ 38 of radio-quiet 7-ray pulsars for GLAST with the ratio of 
radio-quiet to radio-loud increasing from 1.5 for EGRET to 3.2 for GLAST. Curvature 
radiation of 7-rays below the curvature pair death line occurs along all open field lines all 
the way out to the light cylinder, as the weak production of pairs through nonresonant inverse 
Compton scattering is insufficient to screen the electric field. The increased sensitivity of 
GLAST results in a large number of detections at higher altitudes where the 7-ray beam is 
at large angles relative to the radio beam, resulting in radio fluxes below the flux thresholds 
of our select group of radio surveys. However, in many cases, the radio fluxes are relatively 
large but are outside of the sky region of the Parkes Multibeam and Swinburne Intermediate 
Latitude surveys. MSPs are much less confined to the plane of the Galactic disk than are 
NPs and therefore often appear at high latitudes, far above the narrow belt covered by these 
more sensitive surveys. Our simulations suggest that GLAST will detect many MSPs as 
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point sources, and only a handful as r adio-loud pulsars. The point s ource s can be identified 



as pulsar candidates by their spectra. iHarding. Usov. &: Muslimovl (120051 ) indicate that the 



curvature radiation energy spectra of MSPs are very hard, with cutoffs in the region of 10 
to 50 GeV, well within the GLAST range. Such detected candidates can easily be followed 
up with pointed radio observations to obtain the periods of the pulsars. Thus, GLAST with 
radio follow-up observations offers the potential to discover many new MSPs in the Galaxy. 

We are very grateful to Isabelle Grenier and Jean-Marc Casandjian for providing a 
GLAST-LAT 1 year point source threshold map. We express our appreciation to the anony- 
mous referee for the insightful comments, which led to a much-improved manuscript. We are 
indebted to Dune Lorimer who pointed out our neglect of the Shklovskii effect in an earlier 
version of the manuscript. We express our gratitude for the generous support of the National 
Science Foundation (REU and AST-0307365), the Michigan Space Grant Consortium and 
the NASA Astrophysics Theory Program. 
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Table 1: Detected and Simulated MSPs in each Radio Survey 



Survey 


Frequency (MHz) 


Detected 


Simulated 


Arecibo 3 


430 


4 


8 


Arecibo 2 


430 


2 


1 


Greenbank 3 


390 





3 


Greenbank 2 


390 


1 





Molongo 2 


408 





1 


Parkcs 2 


436 


18 


23 


Parkes 1 


1520 








Jodrell Bank 2 


1400 








Parkes MB 


1374 


28 


20 


Swinburne IL 


1374 


12 


9 
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Table 2: Sim ulated 7-ray MSP statistics 

Instrument Radio- loud Radio-quiet Ratio RQ/RL 
EGRET detected 1 ? 

EGRET simulated 4 5-6 1.5 

AGILE simulated 7 11-13 1.9 

GLAST simulated 12 33 - 40 3.2 
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Fig. 1. — Height above the Galactic plane Z (a) and Galactic radius R (b) distributions 
of millisecond pulsars. Initial distributions are indicated by dotted curves and equilibrium 
distributions by solid curves. 
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Fig. 2. — P — P diagram of radio pulsars (dots) in the ATNF catalogue 
(littp://www.atnf.csiro.au/people/pulsar/psrcat/ and iManchester et al.l (120051 ) ). Millisec- 
ond pulsars detected by the group of ten radio surveys used in this study are indicated with 
an additional open circle. The thin dashed lines represent the indicated traditional magnetic 
surface strength, assuming a constant dipole spin-down field. The dotted broken lines repre- 
sent the pair death lines for curvature radiation (CR) and for nonresonant inverse Compton 
scattering (NRICS). The upper and lower MSPs birth lines for Eddington critical accretion 
rate Me (above) and 5 x 1{)~^Me of that rate (below) discussed in the text are indicated by 
solid lines. 




Fig. 3. — Geometry defining the angular quantities associated with a charged particle along 
an open field line defined by ^ emitting curvature radiation at 6 and r along the line of sight, 
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Fig. 4. — Aitoff projections of detected (a) and simulated (b) MSPs. 
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Fig. 5. — P — P diagram of millisecond pulsars detected (a) and simulated (b) by the group of 
ten radio surveys. The thin dashed lines represent the indicated traditional magnetic surface 
strength, assuming a constant dipole spin-down field. The dotted broken lines represent the 
pair death lines for curvature radiation (CR) and for nonresonant inverse Compton scattering 
(NRICS). The upper and lower MSPs birth hnes for Eddington critical accretion rate Me 
(above) and 5 x 10"^ Me of that rate (below) discussed in the text are indicated by solid 
hnes. 
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Fig. 6. — Distributions of various characteristics indicated for detected (shaded histograms) 
and simulated (unshaded histograms) MSPs from the Galactic disk. Also indicated is the 
p-value of the Kolmogorov-Smirnov test of the binned detected and simulated sample dis- 
tributions at a significance level of a = 5%. The dotted histograms represent the simulated 
distributions without the Shklovskii effect. 
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Fig. 7. — P — P diagram of radio-quiet (crosses) and radio-loud (solid circles) 7-ray normal 
and millisecond pulsars detected by EGRET (upper left) and millisecond pulsars predicted 
to be detected by EGRET (upper right), AGILE (lower left) and GLAST (lower right). The 
thin dashed lines represent the indicated traditional magnetic surface strength, assuming a 
constant dipole spin-down field. The dotted broken lines represent the pair death lines for 
curvature radiation (CR) and for nonresonant inverse Compton scattering (NRICS). The 
upper and lower MSPs birth lines for Eddington critical accretion rate Me (above) and 
5 X 10~^ Me of that rate (below) discussed in the text are indicated by solid lines. 
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Fig. 8. — Radio flux at 1400 MHz plotted for detected and simulated radio pulsars in (a) 
and for GLAST radio-quiet and radio- loud 7-ray pulsars in (b). 
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Fig. 9. — The impact angle plotted as a function of the radio flux at 1400 MHz for the 
simulated GLAST radio-quiet 7-ray pulsars (10 times). 
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Fig. 10.— Log(N)-Log(S) for the radio flux at 400 MHz (a), 1400 MHz (b) and for EGRET 
(c) and GLAST (d). The hght histograms represent the simulated undetected source counts. 
For the radio flux, the dark histograms indicated the simulated detected flux at 400 MHz (a) 
and 1400 MHz (b). For EGRET (c) and GLAST (d), the predicted 7-ray source counts are 
indicated for radio-loud (dark histograms) and radio-quiet (medium dark histograms) 7-ray 
millisecond pulsars. 



